/*
 *  SurfaceEllipsoid.cpp
 *  EPPI-0.0
 *
 *  Created by Roberto Ortega on 18/04/11
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "surfaceellipsoid.h"

// ***
SurfaceEllipsoid::SurfaceEllipsoid()
:Surface("meshname")
{}

// ***
SurfaceEllipsoid::~SurfaceEllipsoid()
{}

// ***
SurfaceEllipsoid::SurfaceEllipsoid(ifstream &in, const string& na)
:Surface(na)
{
	rx = 1.0;
	ry = 1.0;
	rz = 1.0;
	
	nn = 10; 
	mm = 10;
	
	string keyword;

	while (in >> keyword && keyword.compare("ENDSURFACE") )
	{		
		if (!keyword.compare("DIMS"))
		{
			in >> rx;
			in >> ry;
			in >> rz;
			
			cout << keyword << "\t";			
			cout << rx << "\t" << ry << "\t" << rz << endl;
		}
		if (!keyword.compare("GRID"))
		{
			in >> nn;
			in >> mm;
			
			cout << keyword << "\t";
			cout << nn << "\t" << mm << endl;
		}		
		if (!keyword.compare("CONTACTPOINTS"))
		{
			in >> pointsOption;
			
			if (!pointsOption.compare("all"))
			{
				contactAllOption = true;
				
				cout << keyword << "\t";
				cout << pointsOption << endl;
			}
			
			if (!pointsOption.compare("file"))
			{
				in >> contactFileName;
				contactFileOption = true;
				
				cout << keyword << "\t";
				cout << pointsOption << "\t" << contactFileName << endl; 
			}
		}		
	}
	
	createEllipsoid();
	this->createContactPoints();	
}

// ***
void SurfaceEllipsoid::createEllipsoid(void)
{
	int i;
	int j;
	
	double x;
	double y;
	double z;
	
	double u;
    double v;
	
	int elem;
    int node;
	
	vector<double> vCoordinates;
    vector<int> vElement;
	
	vCoordinates.resize(3);
	
	for (j = 0; j < mm + 1; ++j)
	{
		for (i = 0; i < nn + 1; ++i)
		{
			u = (pii * i) / nn;
			v = (2 * pii * j) / mm;
			
			node = i + j * (nn + 1);
			
			x = rx * cos(v) * sin(u);
			y = ry * sin(v) * sin(u);
			z = rz * cos(u);
			
			vCoordinates[0] = x;
			vCoordinates[1] = y;
			vCoordinates[2] = z;
			
			vLocalPoints.push_back(vCoordinates);
		}
	}
	
	vElement.resize(4);
	
	for (j = 0; j < mm; ++j)
	{
		for (i = 0; i < nn; ++i)
		{
			elem = i + j * nn;
			
			vElement[0] = elem + j;
			vElement[1] = elem + j + 1;
			vElement[2] = elem + j + nn + 2;
			vElement[3] = elem + j + nn + 1;
			
			vLocalElements.push_back(vElement);
		}
	}
	
	nLocalPoints = vLocalPoints.size();
	nLocalElements = vLocalElements.size();
}

// ***
double SurfaceEllipsoid::evalS(double& x, double& y, double& z)
{
	S = pow(x / rx, 2) + pow(y / ry, 2) + pow(z / rz, 2) - 1.0;
	
	return S;
}

//***
Eigen::RowVector3d SurfaceEllipsoid::evaldS(double& x, double& y, double& z)
{
	dS.x() = 2.0 * x / (rx * rx);
	dS.y() = 2.0 * y / (ry * ry);
	dS.z() = 2.0 * z / (rz * rz);
	
	return dS;
}

//***
Eigen::Matrix3d SurfaceEllipsoid::evalddS(double& x, double& y, double& z)
{	
	ddS(0,0) = 2.0 / (rx * rx);
	ddS(1,1) = 2.0 / (ry * ry);
	ddS(2,2) = 2.0 / (rz * rz);
	
	return ddS;
}